% To calculate serial dependence of up- and downgrade for each Pij of std
% transition prob (cohort method) of G1000k annual data
% Reference: Gagliardini and Gourieoux (2005)



clear all;
clc;
load Prob_std_TP06YRJfinR5G1000kMR7.mat

TM =zeros(12,7,8); % take the needed data from prob_std. totaly 12yrs from 1998-2009. 7 ratings and 8 columns including default
TM(:,:,1:7)=prob_std(2:13,3:9,3:9);
TM(:,:,8)=prob_std(2:13,3:9,11);

TimeLength = size(TM,1);
RateSize = size(TM,2);


% up- and downgrade prob involving a migration of at most two buckets
% downgrade porb: d(k,t) = p(k,k+1,t) + p(k, k+2,t)
% upgrade porb:   u(k,t) = p(k,k-1,t) + p(k, k-2,t)
% except for the extreme categories, where the number of buckets is one)

DownProb = zeros(TimeLength, RateSize); 
for yr = 1: TimeLength
    for rate = 1: RateSize % last grade 4.5 can downgrade to default
        if rate == RateSize
            DownProb(yr,rate) = TM(yr, rate, rate+1);
        else
            DownProb(yr,rate) = TM(yr, rate, rate+1)+ TM(yr, rate, rate+2);
        end
    end
end

UpProb = zeros(TimeLength, RateSize-1);
for yr = 1: TimeLength
    for rate = 2: RateSize
        if rate == 2
            UpProb(yr,rate) = TM(yr, rate, rate-1);
        else
            UpProb(yr,rate) = TM(yr, rate, rate-1)+ TM(yr, rate, rate-2);
        end
    end
end
 
% figure for up- and downgrade prob along time
year_list=[1998:2009]';
year_list_text=num2str(year_list);
Y = cellstr(char('1998','2000','2002','2004','2006','2008'));
% downgrade
figure;
for i = 1:RateSize
subplot(3,3,i); plot (DownProb(:,i),'-o');
set(gca,'XTick',1:2:TimeLength);
set(gca,'XTickLabel',Y);
RateText = strcat('Rating', num2str(1+0.5*i));
legend(RateText,'Location','NorthEast');
end

% upgrade
figure;
for i = 2:RateSize
subplot(2,3,i-1); plot (UpProb(:,i),'-s');
set(gca,'XTick',1:2:TimeLength);
set(gca,'XTickLabel',Y);
RateText = strcat('Rating', num2str(1+0.5*i));
legend(RateText,'Location','NorthEast');
end

% autocorrelation order 1 and 2
% http://www.itl.nist.gov/div898/handbook/eda/section3/eda35c.htm

% using corr(X,Y) return 7x7 matrix which contains the corrs between each rating
% we only need the diag elements
DownCof_order_1 = diag (corr(DownProb(1:end-1,:), DownProb(2:end,:))); 
DownCof_order_2 = diag (corr(DownProb(1:end-2,:), DownProb(3:end,:)));
UpCof_order_1 = diag (corr(UpProb(1:end-1,:), UpProb(2:end,:)));
UpCof_order_2 = diag (corr(UpProb(1:end-2,:), UpProb(3:end,:)));


